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In this paper, we extend our analysis of lattice systems using the wavelet transform to systems 
for which exact enumeration is impractical. For such systems, we illustrate a wavelet-accelerated 
Monte Carlo (WAMC) algorithm, which hierarchically coarse-grains a lattice model by computing 
the probability distribution for successively larger block spins. We demonstrate that although the 
method perturbs the system by changing its Hamiltonian and by allowing block spins to take on 
values not permitted for individual spins, the results obtained agree with the analytical results 
in the preceding paper, and "converge" to exact results obtained in the absence of coarse-graining. 
Additionally, we show that the decorrelation time for the WAMC is no worse than that of Metropolis 
Monte Carlo (MMC), and that scaling laws can be constructed from data performed in several short 
simulations to estimate the results that would be obtained from the original simulation. Although 
the algorithm is not asymptotically faster than traditional MMC, because of its hierarchical design, 
the new algorithm executes several orders of magnitude faster than a full simulation of the original 
problem. Consequently, the new method allows for rapid analysis of a phase diagram, allowing 
computational time to be focused on regions near phase transitions. 



I. INTRODUCTION 



One of the fundamental challenges in the simulation of very large systems is balancing the competing aims of 
accuracy, both numerical and physical, and computational efficiency. As the number of degrees of freedom, the 
number of interactions, and the number of system parameters become large, the computational cost and storage 
requirements of any algorithm which models these systems rapidly become prohibitive. For a system with TV degrees 
of freedom, the complexity of simulation algorithms is typically O {N'^) , although this can be reduced to O {N) for 
methods incorporating both cell and Verlet lists^ and O (iV"^/^) for methods including Ewald sums,^ or increased to 
O {N^^ or more for quantum methods.^ Consequently, very large systems are still expensive to simulate, even with 
efficient algorithms. One way to reduce the complexity of such systems is to employ a coarse- graining method, which 
systematically reduces the number of degrees of freedom in the system, and thereby the overall complexity of the 
simulation. Coarse-graining techniques have been developed for both on- and off-lattice calculations. Lattice simula- 
tions have generally relied on coarse-grainings based on renormalization group theory,*"^ while off-lattice simulations 
include coarse-graining techniques as varied as united-atom models, mapping to lattice models,^' ^ dissipative particle 
dynamics,^ and bead-and-spring models. 

At the same time, any benefits obtained from applying a suitable coarse-graining technique must be weighed against 
the principal difficulty of using such a method — inaccurate numerical results, such as a different critical point. The 
use of uncontrolled approximations also presents a problem, since even if a simulation converges to give an answer, we 
often cannot use results from those simulations to provide any insight into the systems about which we are actually 
interested. In addition, in virtually all cases the transformation is irreversible: we cannot reconstruct our original 
system after we have "evolved" a coarse-grained system. 

The work presented in this paper represents the application of a new technique to enhance the performance of 
traditional lattice simulations using the wavelet transform method. The theoretical foundations of this technique 
are outlined in the preceding paper (henceforth denoted as I^^). First developed as an analysis technique which has 
found wide acceptance in signal processing, the wavelet transform has not yet been extensively applied to the field of 
molecular simulations. Exceptions to this include the variational work of Best and Schafer,^^' ^'^ applied to statistical 
field theory, as well as the work of Huang in studying self-similarity in high-energy physics. 

We choose to study the class of Ising lattices with pairwise and external field interactions, so that the Hamiltonian 
can be written in the form 

W = -^/ijO-i - JifcCTiO-fc, (1) 

i i k 

where the indices i and k run over all spins in the system, and the interaction strengths hi and J,;fc can vary with 
position in the lattice. The nearest-neighbor Ising lattice corresponds to Jik = unless Oi and Ufc are neighboring 
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spins on the lattice, in which case Jik = J. Our goal is to construct a phase diagram for such systems as efficiently 
as possible. The original problem would be studied using the Metropolis Monte Carlo (MMC) technique, along 
with possible improvements such as the Swendsen-Wang or Wolff spin-cluster algorithms. We introduce a new 
method for studying this problem, using the wavelet transform technique to simulate the system hierarchically, by 
systematically transforming (1) to provide expressions for potentials at coarser scales. 

Since the systems under consideration in this paper are lattices, we present the wavelet transform using the discrete 
framework, involving filter banks and matrices. In the discussion below, we present equations for one-dimensional 
systems only; however, the simulations are performed using higher-dimensional wavelets created using the lifting 
scheme of Swcldcns.^** The lifting scheme can be applied recursively as needed to achieve representations of the data 
at multiple length scales. However, we do not move between different physical models for the system, but only change 
the length scales which the different variables describe. As a result, the wavelet transform can be described as a 
multiresolution tcchnique;^*^^^ multiresolution techniques comprise a subset of the more general class of multiscale 
techniques, which can describe any simulation involving multiple modeling steps, possibly involving the use of multiple 
underlying physical models. We refer the reader to paper I for further details on the implementation of the wavelet 
transform with respect to lattice systems. 



II. THE WAVELET-ACCELERATED MONTE CARLO (WAMC) ALGORITHM 

The principal difficulty of performing a wavelet transformation on a lattice system is in working with the discrete 
set of values that each spin is permitted to take, such as in a spin-g Ising model. Because the transformed variables 
are linear combinations of the original variables, the constraint that the spins of the individual sites on the original 
lattice must be drawn from the set {— (j, —q + 1, . . . , (j — 1, (7} quickly becomes a more complicated constraint on the 
transformed variables Sj and 5j. As the system size becomes large, the difficulty of rewriting the spin constraints 
proves so great that previous investigations of the use of wavelets in statistical field theory ignored Ising models 
altogether. '^"^ Consequently, we would like, if at all possible, to avoid computations involving original states after we 
have carried out the wavelet transformation. At this point, we take note of the application of wavelets to image 
compression, where the goal is to reduce the amount of information needed to reconstruct an image. We would like 
to apply this technique to lattice systems, and reduce the number of degrees of freedom which must be accounted for 
in our calculations. 

We consider our system to be a d-dimensional regular lattice C with side length I, so that the size of the lattice is 
= |£| = Z"^, and we let a site CTj on the lattice C be characterized by a "spin" chosen from a finite set J of values and 
by its physical location on the lattice. For the spin-i Ising model, for example, the set J is just { + |, — ^} (although 
for computational convenience this is usually treated as {+1,-1}, a convention which we follow below); similarly, for 
a lattice gas based on a spin-1 Ising model, J = {0, 1,2} represents the allowed occupation numbers of each lattice 
site. We then assume that the only physical interactions that occur are either interactions with an external field hi 
which can vary at each lattice site, or pairwise interactions with the bilinear form U ((Ti,o"j) = JijCJiGj, where Jij is 
usually a function only of the spacing between sites i and j. Consequently, the Hamiltonian of the system can be 
written in the form (1), 

i i j 

For the purposes of our simulations, however, we will find it more convenient to treat the set of spins (ci, . . . ,<TAr) 
and the external field {hi, . . . , hpf) as vectors u and h, and the pairwise interaction strengths Jij as a matrix J. Then 
the Hamiltonian (2) can be written in matrix form as 

-f3n = h^u + u^Ju. (3) 

This formulation of the problem is similar in spirit to that of graph theory, where the pairwise potential Jij is used to 
generate an adjacency list which specifies which edges interact. Using (3) as the basis for a Monte Carlo simulation 
requires the calculation of the change of energy AEnm from microstate to microstate u„: 

AEnm = h"^ (U„ - Um) + (u„ - U^)"^ J (u„ - U„) . (4) 

If moves are restricted to changes of single spin flips, then only a single entry of u„ — is nonzero, and therefore 
the calculation (4) reduces to a dot product, instead of a matrix multiplication. 

As described in I, the action of the wavelet transform is to insert between each product in (3) or (4) the identity 
matrix in the form I = W^W, where W is the wavelet transform which maps data from one scale to the next 
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coarser scale, containing half as many data points. The resulting expressions rewrite the Hamiltonian in terms of 
wavelet-transformed averages and differences, with downsampling needed to reduce the number of variables from 2N 
to A''. As before, the wavelet transform can be iterated by applying it to successive sets of averages, leading after K 
iterations to Hamiltonians of the form: 



where in (5) the u*^^^ represent "block spins" whose values are determined by wavelet averaging over some well-defined 
region of the original system. The Hamiltonians (3) and (5) have the same formal structure, so that Monte Carlo 

simulations of the two systems are essentially identical. The only modifications needed to simulate a coarse-grained 

Hamiltonian are the ability to select new microstates u^^' which are generated through wavelet transformations of 
the original microstates Uj, and the elimination of unwanted degrees of freedom from (5). It should be noted that in 
(5), the elements of ii^^'' are not restricted to the same values as in the original system, but are free to take on any 
value which is consistent with the wavelet transform applied to the system. 

Prom above, we saw that for even Hamiltonians H (x), we should have that (^)^ = for any wavelet difference 5, 
where (•) ^ denotes the ensemble average weighted by the Hamiltonian H. As a "worst-case scenario" for our method, 
we shall assume not only that {5) jj = 0, but also that any terms in the Hamiltonian (5) containing fluctuation terms 
can be neglected as well. This assumption allows us to reduce the size of J(^) from NxN to x where 

d is the lattice dimensionality. Consequently, instead of performing calculations involving all of the original variables x 
which describe the state of our system, we consider functions only of local averages of our original variables. However, 
we anticipate that this simplification of the interactions present in the system will have a significant impact on the 
thermodynamic behavior of the resulting system; we illustrate these effects below. 



To generate the new microstates u^- , we need an estimate for the probability distribution p [u^ ] which describes 



the individual sites in the coarse-grained lattice. Determining the correct distribution for a given ' would require 
a detailed simulation of the original system. An alternative, ignoring the effect of neighboring block spins, would 
be to perform an exact enumeration of the spins within a block, which is possible only for the smallest of block 
spins. Since we would like to apply this method to systems of arbitrary size, we want to avoid both of these options. 
Therefore, we simulate a sublattice with the same dimensions as u , ignoring physical interactions with the rest of 
the system by using either free or periodic boundary conditions. Using the standard Metropolis acceptance criterion, 
we compute distributions for the properties of the small lattice, such as the magnetization. Then, according to the 
matrix formulation described in Sections H and HI of Paper I,^^ since the wavelet transform defines a single block spin 
ui^^ as a linear function of the individual spins at level K — 1 which it replaces, we can use the linearity properties 
of probability distributions to convert the distribution of the properties directly into a distribution for the block spin 
^(^) 25 Finally, using the distribution for the block spin iif^^ as a starting point, we perform a Monte Carlo simulation 
on the system of block spins defined by the Hamiltonian (5). 

Although (3) and (5) are structurally the same, we cannot impose a one-to-one correspondence between the states 
in the configuration space of (4) and the states in the configuration space of (5). Consequently, the thermodynamic 
information obtained from the two will not necessarily be identical; as we have shown in paper I, there is under fairly 
broad conditions a loss of entropy associated with the application of coarse-graining to a system. We can ensure that 
the detailed balance condition for the simulation based on (5) is satisfied for the new simulation by requiring 



where a (m — > n) is the probability of accepting a move from microstate m to microstate n, and p (m) is the probability 
of selecting microstate m as determined from simulations on finer-grained lattices at lower scales. 



The wavelet transform is a hierarchical method which can be applied iteratively to a system to obtain successively 
coarser descriptions of a system. To describe the operation of the wavelet transform on a lattice model, we need to 
introduce some notation based on the various length scales in the problem. In the original problem, the applicable 
length scales are the lattice spacing /, the correlation length ^, and the total lattice size L. Applying the wavelet 




(5) 
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transform method once increases the lattice spacing by some factor a, so that the ratios of correlation length to lattice 
spacing and of system length to lattice spacing each decrease by a. If we apply the wavelet transform m times in 
succession, the corresponding factor becomes a™. 

We perform the simulation in a series of K stages, where the length scales at each stage are functions of the length 
scales at the previous stages. The initial simulation is performed on a sublattice of the original problem, with lattice 
size L^-"^^ < L, where the superscript denotes the first stage of the simulation. The lattice spacing of the first stage 
is the same as in the original problem, so we define l^^^ = I. At each subsequent stage of the simulation, the lattice 
spacing of the fcth stage is defined by the recursive relation l^''^ = L^''~^H^''~^\ Since l^^^ is fixed to be the lattice 
spacing of the original lattice, the adjustable parameter in this relation is the lattice size i'^'^^ of each stage. If we 
assume that the lattice is the same length in all dimensions at every stage, a single variable in stage /c is a block 

variable representing the (i^*^"^)) variables simulated in stage k — 1. 

Assimiing that the lattice is the same length in all directions both in the original problem and at every stage in 
the wavelet-transformed problems, there are Nt = L'^ lattice variables in the original problem, and iV^*^' = (L^*^^) 
lattice variables in the k*'^ stage of the wavelet-transformed problem. However, each variable in stage fc is a block 
variable representing the average behavior of the (L^*^^^^) variables in a block at stage fc — 1, so the number of total 
degrees of freedom represented at stage k is N^''^ = HiLi where N^^^ = Nf. The number of simulated degrees 
of freedom is TV,, = Ni for traditional Metropolis Monte Carlo (MMC), but TV,, = E^Li for WAMC. Because 
the running time of Monte Carlo simulations is usually linear in the number of degrees of freedom being simulated, 
the advantage of coarse-graining the system using a wavelet transform becomes evident. For example, consider an 
"original problem" of simulating a cubic lattice with 256 Ising variables on a side. If we divide the original problem 
into two stages consisting of cubes of 16 Ising variables on a side, we reduce the original problem of analyzing 256^ = 
16,777,216 variables to the simpler problem of analyzing 2 (l6'^) = 8192 variables. Although it is more difficult to 
produce a trial configuration in a simulation of the wavelet-transformed problem than in a simulation of the original 
Ising lattice problem, this is more than offset by the reduction in the number of degrees of freedom being simulated. 

IV. RESULTS 

For the purposes of comparison, our "experimental" systems are two-dimensional Ising models of size 32 x 32, where 
we have run both MMC simulations on the full lattice, and WAMC simulations at a variety of resolutions; we shall 
denote these resolutions using the notation {x, y), where x indicates the length of the block size simulated in the first 
stage to estimate the probability distribution p (u^^') to be used in the second stage, and y denotes the number of 
blocks on a side of the lattice in the second stage of the simulation. 

A. Order parameter 

Usually, the property of greatest interest in a simulation of a lattice system is the order parameter q. For spin 
systems, r] is generally taken to be the magnitude of the average magnetization, so that 77 = (m). [For XY and 
Heisenberg models, and other models where spins are oriented, we generally consider only the magnitude of the 
average vector 77 = (m) = (|m|).] Generally, this is a very simple property to compute, since the value of the order 
parameter is constantly updated during the course of the simulation, and is thus always available. 

For the 32 x 32 Ising model, the results of a MMC simulation, as well as (4, 8)- and (8, 4)-WAMC simulations are 
shown as Figure 1. The primary difii'erence in the curves for the three cases is that as the coarse- graining process 
decreases the number of degrees of freedom in the final stage of the simulation, the location of the Curie temperature, 
indicating onset of spontaneous magnetization, increases and the steepness of the curve below the Curie temperature 
decreases. This result is consistent with our findings for average absolute magnetization (|m|) from analytical models, 
discussed in paper I. In the present case, we note further that we achieve agreement between the different models 
not only in the low-temperature region, but also in the high-temperature regime T ^ Tc- The differences in the 
intermediate regime can be attributed largely to the difference in behavior that rcsidts from the use of the wavelet 
transform to move from the original Hamiltonian (3) to a coarse-grained Hamiltonian (5). Additionally, the increased 
noise in the WAMC results at intermediate and high temperatures arises because of the approximations used for the 
probability distributions p (u^^^) at the second stage of the simulation. The relative lack of noise in the MMC results 
stem in part from the fact that the Metropolis technique leads to non-ergodic sampling of phase space as temperature 
increases, as the simulation tends to cycle through a limited number of states. 
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B. Internal Energy 



Plotting the internal energy {U) as a function of the temperature, we obtain curves that follow the same general 
pattern outlined in paper I. As illustrated in Figure 2, at low temperatures, the internal energy, as computed for 
the 32 X 32 model using standard MC as well as (4, 8)- and (8, 4)-WAMC simulations, is in exact agreement for all 
methods. This occurs because only a few microstates of the system, corresponding to states that have all spins aligned, 
arc actually observed by the system, and the wavelet transform preserves the energy of these states exactly. All three 
eventually reach an average internal energy of zero, but exact agreement is only expected in the infinite-temperature 
limit, when the difference in energy levels between microstates becomes unimportant. For intermediate temperatures, 
as before, the disagreement is a result of the change in form of the Hamiltonian that results from neglecting local 
correlations. Also, we note that for WAMC the "noise" in the internal energy increases both with increasing proximity 
to the "observed" critical point of the system as well as with inc;reasing coarse graining. The additional coarse graining 
yields a Hamiltonian with reduced numbers of energy levels, since the energy of a block spin is defined here to be a 
function only of its overall magnetization, and not of its internal magnetization fluctuations; the reduced number of 
discrete energy levels yields noisier data. 

C. Fluctuation properties 

Fluctuation properties are useful for locating critical points since in the vicinity of a critical point, the magnitude 

of fluctuation properties is known to diverge as where t = \T — Tc \ /Tc.^^ Thus, a rapid increase in the value of a 

fluctuation property such as the heat capacity at constant external field, Ch = {^{E"^) ~ {E}"^^ /ksT^, with respect 

to temperature can be used to estimate the critical temperature of a system. However, use of the wavelet transform 
leads to a decrease in the magnitude of the heat capacity, since the coarse-graining leads to smaller variances in the 
distribution of the energy (U). Consequently, the maximum value of the heat capacity Cmax decreases as a function 
of the number of degrees maintained in the problem. 

In Figure 3, the heat capacity is shown as a function of dimensionless temperature ksT/J for the same systems as 
for the order parameter and internal energy measured above. We see that the location of the maximum of the heat 
capacity does increase, as expected. Although the relative maximum of the heat capacity obtained from the (4, 8)- 
and (8,4)-WAMC simulations appears to be identical, they differ by about 3 percent. Moreover, the actual value of 
the maximum is not as important as its existence and its location as a function of ksT/ J and of the resolution of the 
model. 



D. Scaling results 

One application of the wavelet-accelerated MC method is to provide an upper bound for locating phase transitions. 
Running multiple simulations, at different levels of resolution, one can determine for cac;li level of resolution the 
approximate phase transition temperature Tp (Ns), where Ns is the number of degrees of freedom (here, block spins) 
in the given system. From these data, one can extrapolate a scaling relationship of the form 

{Tp - t;) ^ N-\ (6) 

where 7 is the corresponding "scaling" exponent, and T* is the phase transition temperature for the untransformed 
model. The estimate obtained for the scaling exponent 7 depends upon the technique used to calculate the transition 
temperature for a given system for example, estimating divergence of the heat capacity versus the onset of sponta- 
neous magnetization. Our simulations suggest that 7 is typically between 0.20 and 0.25, with lower values obtained 
from divergence of heat capacity than from the onset of spontaneous magnetization. 

As explained in I, using a relationship like (6) to estimate the phase transition temperature will usually lead to an 
overestimate of the phase transition temperature. This is a consequence of the underestimation of entropy that occurs 
through the reduction of the size of configuration space as a result of the wavelet transform. Estimates of ksTp/J, 
as determined by (6) for the two-dimensional Ising model considered here typically varied between about 2.7 and 2.9, 
which is an error of approximately 25 percent from the theoretical value of 2.27 provided by the Onsager solution, but 
only about 20 percent from the results determined by the traditional MC simulations, which gave ksTp/J « 2.35. 
Although these errors are somewhat sizable, it is useful to note that the total computation time required to obtain the 
estimate using scaling laws is at least an order of magnitude smaller than the computation time required to perform 
a direct simulation on the original system. Thus, if computational time is at a premium, an effective approach may 
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be to use the wavelet transform method to provide an upper bound for Tp, and then perform a direct simulation for 
the parameter space with temperatures below Tp. 

E. Decorrelation time 

Another important measure to study is the time required for decorrelated samples. It is well known that in the 
vicinity of the critical point, traditional Monte Carlo algorithms experience so-called "critical slowing-down."^'^ The 
Monte Carlo aspect of the WAMC algorithm does not vary from traditional Metropolis Monte Carlo, so we expect that 
the performance of the two algorithms should be similar, when measured near their respective critical temperatures. 

To compare the two methods, we generated 2^^ = 16 777 216 new configurations for the 32 x 32 Ising model at the 
critical point using traditional Metropolis Monte Carlo, as well as for the (4, 8)-, (8, 4)-, and (16, 2)-WAMC models. 
To determine the correlation time, we used the blocking technique of Flyvbjerg and Petersen. The results are shown 
as Figures 4 and 5 for the MMC and (16, 2)-WAMC models, respectively. The salient feature in the graph is the 
onset of a plateau in the value of the variance of the energy ct^; according to the method of Flyvbjerg and Petersen, 
this indicates that configurations separated by a distance of 2^ steps are statistically independent, where x is the 
number of blocking transformations that have been performed. For the MMC model, we find that a; « 13 or a; « 14 
provides a decent estimate; for the (16, 2)-WAMC model, a; = 16 is a good estimate for the index. [For the (4, 8)- and 
(8,4)-niodels (not shown), a; = 15 is a reliable estimate] In each case, this indicates that between 2^^ = 16 384 and 
2^^ = 65 536 steps are required between independent configurations. Thus, we conclude that there is no degradation 
of performance near the critical point of a WAMC simulation, relative to traditional MMC simulations. 

V. ANALYSIS 
A. Measured performance comparison 

In comparing the performance of the standard Monte Carlo algorithm to the wavelet-accelerated Monte Carlo 
algorithm, we performed 5 x 10^ lattice passes on a 32 x 32-lattice on a 733 MHz Pentium II: for the standard 
Monte Carlo algorithm, this meant that, on average, 5 x 10^ attempts were made to flip each spin. For the WAMC 
algorithm, 5 x 10^ attempts were made to flip a spin on a given level. The results are summarized in Table 1. We 
see that the (4,8)- and (8, 4)-simulations, which have the smallest total number of lattice sites (80), perform the 
fastest; however, even the (16, 2)-simulation, which has one-fourth as many variables (260) as the 32 x 32 standard 
Monte Carlo simulation (1024), finishes in less than 8 per cent of the time required for the latter simulation. As the 
system size increases, the computational efficiency achieved by breaking down the system into multiple stages, all of 
relatively equal size, becomes even greater: for a 128 x 128-lattice, the performance gain increases from a factor of 
approximately 25 to a factor of approximately 50, when we compare the (32,4)- and (16,8)-WAMC simulations to 
the standard MC model. However, for the (8, 16)-WAMC model, the complexity of assigning one of 65 possible values 
to each of 256 variables according to the correct probability distribution becomes comparable to that of the original 
problem, so that in fact standard MC runs in roughly a factor of 3 faster than the (8, 16)-model. 

Similar results are also observed for calculating the phase diagram of a 64 x 64-Ising lattice, which is shown in 
Figure 6 as a plot of average magnetization as a function of temperature and external field strength using an (8, 8)- 
modcl, for temperatures between T = 0.5 and T = 5.0, and for field strengths between h = —1 and h = 1. The 
phase diagram reproduces the essential features of the original two-dimensional ferromagnetic Ising lattice, such as 
the phase separation at /i = 0, although the exact shape differs from the results obtained via a standard Metropolis 
Monte Carlo simulations. However, the plot based on WAMC calculations is created approximately 40 times faster 
than would a comparable plot using standard MMC. 

B. Comparison with renormalization group methods 

Our observations also indicate that the accuracy of the wavelet-accelerated Monte Carlo simulations depends on 
the relative proximity to an "attractive fixed point" of the physical model in parameter space. Borrowed from 
renormalization group theory, these attractive fixed points represent the limiting behavior of the system under various 
conditions, such as the zero- and infinite-temperature limits and the limits of zero and infinite external field. As we 
approach these limiting cases, the approximations made in obtaining our wavelet-transformed Hamiltonian become 
increasingly less significant. 
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Combining these observations allows us to design an on-line fine-tuning algorithm for the coarse graining of our 
system: the further away from the critical point of the parameter space, the smaller the number of degrees of freedom 
N^^^ necessary to simulate the system must be. Thus, if we keep track of changes in fluctuation properties such as 
the heat capacity or the magnetic susceptibility as we change the system parameters {h, T), we can get an estimate of 
our relative distance to the critical point. If we are sufficiently far away from the critical point, we can choose either 
to increase the number of stages K that we simulate, or we can look at more degrees of freedom at lower stages by 
increasing N^^\ N^^\ . . . , N^-^~^\ As we approach the critical point, we can either decrease the number of stages K 
or include more degrees of freedom at higher stages by reducing N'-^\ . . . , N^^~^\ 

C. Sources of error 

In general, the source of our errors can be traced to the assumption that local fluctuation terms could be reasonably 
ignored in our coarse-grained Hamiltonian (5). This naive but otherwise useful assumption yields correct thermo- 
dynamic behavior when the overall physics of the system is particularly simple: in the low- and high-temperature 
regimes, for instance, when the number of observed microstates is small or when the differences between observed 
microstates is inconsequential. For more complicated behaviors, as found at intermediate temperatures and above 
all in the vicinity of a critical point, the use of this assumption has a drastic effect on both the phase space of the 
system, which in turn affects all thermodynamic properties of the system, including the internal energy and entropy 
of the system, as well as fluctuation properties of the system. 

More complicated methods for dealing with fluctuation terms have significant drawbacks associated with them. 
Treating the fluctuation terms just like the block averages maintained in u^^^ means that the wavelet-transformed 
Hamiltonian is no simpler than the original Hamiltonian, which affords few advantages in computational time. Like- 
wise, other approaches, such as parametrizing the probability distribution for the elements of u^^^ using a property 
like the energy E, introduce new functional dependencies which cannot be taken into account using the wavelet trans- 
form. Thus, we sacrifice one of the major advantages of the method — moving from one level to another is achieved 
exclusively through use of the wavelet transform. Thus, the most promising avenue for dealing with fluctuation terms 
is to develop a probability distribution for the fluctuation terms via the same approach used to determine probability 
distributions for the local averages. Then, using the probability distribution for the fluctuation terms, we can treat 
the discarded terms of the Hamiltonian as a noise term which can be used to restore some of the entropy that was 
lost as a result of the coarse-graining [see paper I for more details]. However, we have presented our results here to 
show what can be achieved under "worst-case" conditions, without the use of inverse coarse-graining methods. 

Many coarse-graining techniques control errors by fltting the parameters of a new Hamiltonian to ensure agree- 
ment with some known structural information about the system, such as the radial distribution function.^"' For 
lattice systems, this iterative approach is reflected in renormalization group theory, and notably the Wilson recursion 
method,^^ which finds the fixed points of the system. As formulated, WAMC creates a coarse-grained Hamiltonian 
by truncating the Hamiltonian obtained after application of the wavelet transform. As an improvement to this, it 
should be possible to use the wavelet transform to determine which terms will appear in the Hamiltonian, and then 
determine the appropriate parameters to ensure the best fit for some desired property of the system using an iterative 
approach. ^^"^^ 

D. Constructing an adaptive algorithm for MC using the wavelet transform 

As pointed out above, the wavelet transform method tends to produce overestimates for the critical point of the 
system; therefore, if we start with the high-temperature limit of our algorithm and slowly reduce the temperature in 
our simulation, we can observe the movement toward the critical point by watching various fluctuation parameters, 

such as the heat capacity Ch = (^{E'^) — i^)^^ /ksT'^- Near the critical point, we expect to see a rapid increase in 

the value of Ch, consistent with the expected logarithmic divergence observed in the limit of flnite-size systems.^' 
If we use the onset of this logarithmic divergence as an indicator, we can then "step down" and use a finer lattice 
including more degrees of freedom. This system will naturally better reflect the physics of our system, particularly 
in the vicinity of the critical point. We expect that very near the critical point, we will have to simulate the system 
at the original scale, since this will be effectively the only level which accurately represents the underlying behavior 
of the system. However, the region of parameter space where this is necessary is relatively small compared to the 
complete parameter space. This is especially true when we consider that as we proceed below the critical temperature 
of the system, the logarithmic divergence of Ch will also vanish. As a result, as we move increasingly far away from 
the critical point, we begin to approach the other fixed-point behaviors associated with the low-temperature limits 
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of the system. Since these arc reasonably well-preserved using the wavelet transformation, we can safely return to 
increasingly coarse-grained descriptions of our system as the simulation proceeds past the critical point. 

As an example, we compute the spontaneous magnetization curve for a 64 x 64 Ising lattice in the temperature 
range 0.5 < T < 10.0, with AT = —0.05, and choosing as our refinement criterion AC^/AT < —0.5, until we 
reach the finest scale, corresponding to the original problem. We begin by coarse-graining the system to an (8, 8)- 
modcl, where we find that the criterion is triggered only at T = 5.1; we then continue with a (4, 16)-modcl, down to 
T = 4.0, at which point the refinement criterion is exceeded. Refining once more, we proceed with a (2,32)-model 
until T = 3.4, at which point the threshold is again crossed. Since the next refinement is the original problem, we 
proceed at this level of resolution until we have passed the critical point, so that AC/z/AT is positive. As a coarsening 
criterion, we select for simplicity the opposite of the refinement criterion, ACjj/AT < 0.5. Using this criterion, we 
find that we coarsen the model to the (2, 32)-, (4, 16)- and (8, 8)-models at temperatures of T = 1.75, T = 1.65, and 
T = 1.55, respectively. The rapid coarsening of the model results from the higher estimates of the critical point in 
the coarsened models. Since we are well past the critical point, we expect changes in the heat capacity as a function 
of temperature to be relatively small, and thus it is possible to obtain accurate results from a relatively coarse model. 
Computationally, the time required to create this diagram was only 28 per cent that required to perform a standard 
Metropolis Monte Carlo simulation with the same number of steps. Moreover, in the regions that were not simulated 
using MMC; the computation time required was just 8 per cent of the time required for MMC. The resulting plot 
of magnetization versus temperature, shown as Figure 7, compares favorably to the analytical solution of Onsager, 
which is also shown. 



VI. CONCLUSIONS 

The WAMC algorithm can dramatically reduce the time required to calculate the thermodynamic behavior of a 
lattice system; the trade-off for these savings in time is in the accuracy of the results obtained, a general feature of 
coarse-graining techniques. The error that results is a function of position in parameter space: the results obtained 
are generally accurate in the vicinity of fixed attractors of the system, and decrease as one approaches critical points 
of the parameter space. Near critical points, deviations from results performed on the original lattice system are the 
result of coarse-graining the Hamiltonian by eliminating local fluctuation terms. Consequently, this suggests that a 
hierarchical simulation which uses fluctuation properties such as the heat capacity Ch to gauge proximity to critical 
points in "real time" would yield dramatic savings in the computation time of the behavior of a lattice system over 
a wide region of phase space, as regions of space close to fixed attractors would be simulated at a very coarse scale, 
with full-detailed simulations reserved only for regions of parameter space close to critical points. 
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FIG. 1. Absolute average magnetization as a function of the dimensionless temperature ksT/J for the 32 x 32 Ising model 
computed using standard MG (left curve), a (4, 8)-WAMG simulation (center), and a (8, 4)-WAMG simulation (right). 



FIG. 2. Internal energy as a function of the dimensionless temperature ksT/J for the 32 x 32 Ising model computed using 
standard MG (left curve), a (4, 8)-WAMG simulation (center), and a (8,4)-WAMG simulation (right). 



FIG. 3. Heat capacity as a function of the dimensionless temperature fesT/J for the 32 x 32 Ising model computed using 
standard MG (left curve), a (4, 8)-WAMG simulation (center), and a (8, 4)-WAMG simulation (right). 

FIG. 4. Graph showing variance in the estimate of energy as calculated using the method of Flyvbjerg and Petersen^® for a 
32 x 32 Ising model measured at its critical temperature. 



FIG. 5. Graph showing variance in the estimate of energy as calculated using the method of Flyvbjerg and Petersen for a 
32 X 32 Ising model in a (16, 2)-WAMG measured at the critical temperature determined from the simulation. 

FIG. 6. Phase diagram plotting average magnetization versus temperature and external field strength for a 64 x 64 ferro- 
magnetic Ising lattice, computed using an (8, 8)-model via WAMG. The general features correspond to those that would be 
produced with standard MMG, but require less than 3 per cent of the computational time. 
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FIG. 7. Phase diagram plotting average magnetization versus temperature for a 64 x 64 ferromagnetic Ising lattice, created 
using an adaptive WAMC algorithm, with refinement and coarsening criterion established using the change in heat capacity 

with respect to temperature ACh / AT. The squares represent the simulation results, while the line reproduces Onsager's 
analytical result for the two-dimensional Ising model with zero external field. The simulation used in a given temperature 
region is shown on the plot. 



10 



Table 1. Performance comparison for Metropolis Monte Carlo (MMC) versus Wavelet- Accelerated Monte Carlo 



(WAMC) algorithms 



Simulation 


Time for 5 x lO'' passes (s) 


32 X 32 MMC 


1824.22 


(4, 8)-WAMC 


64.5781 


(8,4)-WAMC 


55.1094 


(16, 2)-WAMC 


144.516 
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